# Función para instalar y cargar paquetesload_packages <-function(packages) {for (pkg in packages) {if (!require(pkg, character.only =TRUE, quietly =TRUE)) {install.packages(pkg, dependencies =TRUE, repos ="https://cran.r-project.org")library(pkg, character.only =TRUE) } }}# Lista de paquetes necesariosrequired_packages <-c("tidyverse", "igraph", "tidygraph", "ggraph", "kableExtra", "DT", "plotly", "visNetwork", "corrplot", "RColorBrewer", "patchwork", "dendextend")# Cargar todos los paquetesload_packages(required_packages)# Configuración globaloptions(scipen =999)theme_set(theme_minimal())# ===============================================# FUNCIÓN DE SIMULACIÓN MEJORADA BASADA EN DATOS REALES CHILENOS# ===============================================simular_datos_chile <-function(n_pacientes =223055, seed =42) {set.seed(seed)# ===============================================# DISTRIBUCIONES EXACTAS DE LOS DATOS REALES# ===============================================# Distribución exacta de sustancia principal dist_principal <-c("Alcohol"=0.3507386071,"Pasta Base"=0.3773329448,"Cocaína"=0.1886216404,"Marihuana"=0.0618546995,"Sedantes"=0.0116832171,"Opioides"=0.0044652664,"Otros"=0.0012104638,"Inhalables"=0.0008069759,"Hipnóticos"=0.0006500639,"Estimulantes"=0.0006007487,"Anfetaminas"=0.0005872991,"Éxtasis/MDMA"=0.0003586559,"Alucinógenos"=0.0002645088,"Crack"=0.0002465760,"Metanfetaminas y Otros Derivados"=0.0001613952,"Tranquilizantes"=0.0001613952,"Lsd"=0.0001255296,"Heroína"=0.0000806976,"Fenilciclidina"=0.0000224160,"Metadona"=0.0000224160,"Esteroides Anabólicos"=0.0000044832 )# Probabilidad exacta de policonsumo por sustancia principal prob_policonsumo <-c("Pasta Base"=0.8439274767,"Alcohol"=0.4834854411,"Cocaína"=0.8360706391,"Marihuana"=0.7993041966,"Sedantes"=0.6266308519,"Opioides"=0.6817269076,"Otros"=0.6333333333,"Inhalables"=0.8333333333,"Hipnóticos"=0.7103448276,"Estimulantes"=0.6791044776,"Anfetaminas"=0.8396946565,"Éxtasis/MDMA"=0.9875000000,"Alucinógenos"=0.9491525424,"Crack"=0.7454545455,"Metanfetaminas y Otros Derivados"=0.9722222222,"Tranquilizantes"=0.7777777778,"Lsd"=0.9285714286,"Heroína"=0.9444444444,"Fenilciclidina"=0.6000000000,"Metadona"=0.2000000000,"Esteroides Anabólicos"=1.0000000000 )# Distribución exacta del número de sustancias dist_n_sustancias <-c("1"=0.2903857793,"2"=0.3405438121,"3"=0.2488579050,"4"=0.1202125036 )# Probabilidades condicionales completas prob_condicionales <-list("Pasta Base"=list("Alcohol"=0.6517239741,"Marihuana"=0.5286576527,"Cocaína"=0.3013687237,"Sedantes"=0.0226694865,"Otros"=0.0136397120,"Inhalables"=0.0097664140,"Anfetaminas"=0.0052871706,"Opioides"=0.0022693249,"Estimulantes"=0.0021980372,"Crack"=0.0016633795,"Alucinógenos"=0.0013307036,"Lsd"=0.0011406031,"Éxtasis/MDMA"=0.0008316898,"Metanfetaminas y Otros Derivados"=0.0008316898,"Hipnóticos"=0.0004990139,"Heroína"=0.0004039636,"Tranquilizantes"=0.0003326759,"Hongos"=0.0002613882,"Metadona"=0.0000831690,"Fenilciclidina"=0.0000237626,"Esteroides Anabólicos"=0.0000118813 ),"Alcohol"=list("Marihuana"=0.2604622031,"Cocaína"=0.2284684408,"Pasta Base"=0.1309788583,"Sedantes"=0.0465015211,"Otros"=0.0214101286,"Anfetaminas"=0.0052406882,"Inhalables"=0.0051000844,"Estimulantes"=0.0045632334,"Opioides"=0.0022880078,"Éxtasis/MDMA"=0.0016872460,"Lsd"=0.0015977708,"Alucinógenos"=0.0011631771,"Hipnóticos"=0.0011631771,"Metanfetaminas y Otros Derivados"=0.0009714446,"Tranquilizantes"=0.0006263261,"Hongos"=0.0005879797,"Crack"=0.0005496332,"Heroína"=0.0003067720,"Esteroides Anabólicos"=0.0001406038,"Metadona"=0.0000511287 ),"Cocaína"=list("Alcohol"=0.7019941530,"Marihuana"=0.4234069356,"Pasta Base"=0.1497635063,"Sedantes"=0.0355097093,"Otros"=0.0151878877,"Anfetaminas"=0.0101252585,"Éxtasis/MDMA"=0.0068927816,"Lsd"=0.0050626292,"Estimulantes"=0.0040405961,"Opioides"=0.0033275497,"Alucinógenos"=0.0029472583,"Inhalables"=0.0023292848,"Metanfetaminas y Otros Derivados"=0.0018539206,"Hipnóticos"=0.0017350795,"Hongos"=0.0011408742,"Tranquilizantes"=0.0010458013,"Crack"=0.0008794239,"Heroína"=0.0004278278,"Esteroides Anabólicos"=0.0001426093,"Metadona"=0.0000713046,"Fenilciclidina"=0.0000475364 ),"Marihuana"=list("Alcohol"=0.6352105530,"Cocaína"=0.3405812858,"Pasta Base"=0.2098282235,"Sedantes"=0.0484163224,"Otros"=0.0232659274,"Lsd"=0.0130463144,"Éxtasis/MDMA"=0.0119591215,"Inhalables"=0.0091324201,"Anfetaminas"=0.0081901863,"Alucinógenos"=0.0071029934,"Estimulantes"=0.0068130753,"Hongos"=0.0050735667,"Opioides"=0.0047836486,"Metanfetaminas y Otros Derivados"=0.0015220700,"Tranquilizantes"=0.0009422338,"Hipnóticos"=0.0007972748,"Crack"=0.0007247952,"Fenilciclidina"=0.0002899181,"Heroína"=0.0002174386,"Metadona"=0.0000724795 ),"Sedantes"=list("Alcohol"=0.4712202609,"Marihuana"=0.2256331543,"Cocaína"=0.1435149655,"Pasta Base"=0.0629316961,"Otros"=0.0376055257,"Opioides"=0.0226400614,"Estimulantes"=0.0072908672,"Éxtasis/MDMA"=0.0072908672,"Hipnóticos"=0.0072908672,"Anfetaminas"=0.0042210284,"Inhalables"=0.0042210284,"Metanfetaminas y Otros Derivados"=0.0023023791,"Hongos"=0.0019186493,"Alucinógenos"=0.0007674597,"Esteroides Anabólicos"=0.0007674597,"Metadona"=0.0007674597,"Tranquilizantes"=0.0007674597,"Crack"=0.0003837299,"Lsd"=0.0003837299 ),"Opioides"=list("Marihuana"=0.3584337349,"Alcohol"=0.3463855422,"Cocaína"=0.1987951807,"Sedantes"=0.1506024096,"Pasta Base"=0.0793172691,"Otros"=0.0301204819,"Estimulantes"=0.0100401606,"Hipnóticos"=0.0100401606,"Anfetaminas"=0.0090361446,"Metadona"=0.0090361446,"Tranquilizantes"=0.0080321285,"Éxtasis/MDMA"=0.0070281124,"Metanfetaminas y Otros Derivados"=0.0070281124,"Hongos"=0.0060240964,"Crack"=0.0030120482,"Fenilciclidina"=0.0020080321,"Lsd"=0.0020080321,"Alucinógenos"=0.0010040161,"Inhalables"=0.0010040161 ),"Otros"=list("Alcohol"=0.4111111111,"Marihuana"=0.3074074074,"Cocaína"=0.2148148148,"Pasta Base"=0.1333333333,"Sedantes"=0.0740740741,"Opioides"=0.0370370370,"Inhalables"=0.0296296296,"Hipnóticos"=0.0148148148,"Éxtasis/MDMA"=0.0111111111,"Tranquilizantes"=0.0074074074,"Alucinógenos"=0.0037037037,"Estimulantes"=0.0037037037 ),"Inhalables"=list("Alcohol"=0.6222222222,"Marihuana"=0.4944444444,"Pasta Base"=0.3222222222,"Cocaína"=0.1000000000,"Otros"=0.0555555556,"Sedantes"=0.0388888889,"Alucinógenos"=0.0111111111,"Opioides"=0.0055555556 ),"Hipnóticos"=list("Alcohol"=0.4620689655,"Cocaína"=0.2413793103,"Marihuana"=0.2068965517,"Sedantes"=0.1034482759,"Pasta Base"=0.0758620690,"Otros"=0.0551724138,"Alucinógenos"=0.0068965517,"Anfetaminas"=0.0068965517,"Estimulantes"=0.0068965517,"Tranquilizantes"=0.0068965517 ),"Estimulantes"=list("Alcohol"=0.4104477612,"Marihuana"=0.3582089552,"Cocaína"=0.2089552239,"Sedantes"=0.0895522388,"Pasta Base"=0.0597014925,"Éxtasis/MDMA"=0.0373134328,"Anfetaminas"=0.0223880597,"Otros"=0.0223880597,"Lsd"=0.0149253731,"Opioides"=0.0149253731,"Metanfetaminas y Otros Derivados"=0.0074626866 ),"Anfetaminas"=list("Alcohol"=0.6335877863,"Marihuana"=0.3664122137,"Cocaína"=0.3129770992,"Pasta Base"=0.1221374046,"Sedantes"=0.0763358779,"Estimulantes"=0.0305343511,"Éxtasis/MDMA"=0.0229007634,"Opioides"=0.0152671756,"Alucinógenos"=0.0076335878,"Lsd"=0.0076335878,"Metanfetaminas y Otros Derivados"=0.0076335878,"Tranquilizantes"=0.0076335878 ),"Éxtasis/MDMA"=list("Marihuana"=0.6875000000,"Alcohol"=0.5250000000,"Cocaína"=0.4875000000,"Sedantes"=0.2000000000,"Opioides"=0.1125000000,"Hongos"=0.1000000000,"Alucinógenos"=0.0625000000,"Estimulantes"=0.0375000000,"Lsd"=0.0375000000,"Otros"=0.0250000000,"Pasta Base"=0.0250000000,"Inhalables"=0.0125000000,"Metanfetaminas y Otros Derivados"=0.0125000000 ),"Alucinógenos"=list("Marihuana"=0.6271186441,"Alcohol"=0.5593220339,"Cocaína"=0.3050847458,"Éxtasis/MDMA"=0.1016949153,"Sedantes"=0.1016949153,"Pasta Base"=0.0677966102,"Hongos"=0.0508474576,"Estimulantes"=0.0338983051,"Lsd"=0.0338983051,"Anfetaminas"=0.0169491525,"Opioides"=0.0169491525 ),"Crack"=list("Marihuana"=0.5454545455,"Alcohol"=0.5272727273,"Cocaína"=0.3090909091,"Pasta Base"=0.1090909091,"Anfetaminas"=0.0181818182 ),"Metanfetaminas y Otros Derivados"=list("Alcohol"=0.6111111111,"Marihuana"=0.3888888889,"Cocaína"=0.1944444444,"Anfetaminas"=0.0833333333,"Estimulantes"=0.0833333333,"Hipnóticos"=0.0833333333,"Pasta Base"=0.0555555556,"Sedantes"=0.0555555556 ),"Tranquilizantes"=list("Alcohol"=0.4722222222,"Marihuana"=0.2500000000,"Cocaína"=0.1111111111,"Estimulantes"=0.0833333333,"Opioides"=0.0833333333,"Otros"=0.0555555556,"Hipnóticos"=0.0277777778,"Lsd"=0.0277777778,"Sedantes"=0.0277777778 ),"Lsd"=list("Marihuana"=0.8214285714,"Alcohol"=0.6428571429,"Cocaína"=0.2857142857,"Sedantes"=0.1071428571,"Anfetaminas"=0.0714285714,"Éxtasis/MDMA"=0.0714285714,"Opioides"=0.0357142857 ),"Heroína"=list("Cocaína"=0.6666666667,"Marihuana"=0.3888888889,"Alcohol"=0.2777777778,"Anfetaminas"=0.1666666667,"Sedantes"=0.1666666667,"Crack"=0.1111111111,"Éxtasis/MDMA"=0.1111111111,"Alucinógenos"=0.0555555556,"Lsd"=0.0555555556,"Metanfetaminas y Otros Derivados"=0.0555555556,"Opioides"=0.0555555556,"Pasta Base"=0.0555555556 ),"Fenilciclidina"=list("Alcohol"=0.6000000000,"Cocaína"=0.6000000000,"Heroína"=0.4000000000,"Éxtasis/MDMA"=0.2000000000 ),"Metadona"=list("Sedantes"=0.2000000000 ),"Esteroides Anabólicos"=list("Alcohol"=1.0000000000,"Cocaína"=1.0000000000 ) )# Para sustancias sin probabilidades condicionales específicas dist_general_otras <-c("Alcohol"=0.35,"Marihuana"=0.25,"Cocaína"=0.20,"Pasta Base"=0.10,"Sedantes"=0.05,"Otros"=0.02,"Opioides"=0.015,"Estimulantes"=0.01,"Hipnóticos"=0.005 )# ===============================================# GENERACIÓN DE DATOS# ===============================================# Inicializar dataframe data_sim <-data.frame(HASH_KEY =paste0("SIM_", sprintf("%08d", 1:n_pacientes)),sustancia_principal =NA_character_,otras_sustancias_no1 =NA_character_,otras_sustancias_no2 =NA_character_,otras_sustancias_no3 =NA_character_,stringsAsFactors =FALSE )# Asignar sustancias principales según distribución exacta data_sim$sustancia_principal <-sample(names(dist_principal),size = n_pacientes,prob = dist_principal,replace =TRUE)# Para cada paciente, determinar policonsumofor(i in1:n_pacientes) { sust_prin <- data_sim$sustancia_principal[i]# Obtener probabilidad de policonsumo para esta sustancia prob_poli <-ifelse(sust_prin %in%names(prob_policonsumo), prob_policonsumo[sust_prin],0.70)# Determinar si tiene policonsumoif(runif(1) < prob_poli) {# Tiene policonsumo - determinar cuántas sustancias en total (2, 3 o 4) probs_n <- dist_n_sustancias[2:4] probs_n <-as.numeric(probs_n) /sum(as.numeric(probs_n)) n_total <-sample(2:4, size =1, prob = probs_n) n_adicionales <- n_total -1# Obtener probabilidades condicionales para esta sustancia principalif(sust_prin %in%names(prob_condicionales)) { probs <- prob_condicionales[[sust_prin]] } else { probs <- dist_general_otras }# Verificar que hay sustancias disponiblesif(length(probs) >0) {# Convertir a vector con nombres probs_vec <-unlist(probs)# Eliminar la sustancia principal si aparece en las opciones probs_vec <- probs_vec[names(probs_vec) != sust_prin]# Eliminar sustancias con probabilidad 0 o NA probs_vec <- probs_vec[!is.na(probs_vec) & probs_vec >0]if(length(probs_vec) >= n_adicionales) {# Seleccionar sustancias adicionales según probabilidades sust_seleccionadas <-sample(names(probs_vec),size = n_adicionales,prob = probs_vec,replace =FALSE)# Asignar a las columnas correspondientesif(length(sust_seleccionadas) >=1) { data_sim$otras_sustancias_no1[i] <- sust_seleccionadas[1] }if(length(sust_seleccionadas) >=2) { data_sim$otras_sustancias_no2[i] <- sust_seleccionadas[2] }if(length(sust_seleccionadas) >=3) { data_sim$otras_sustancias_no3[i] <- sust_seleccionadas[3] } } } } }return(data_sim)}# ===============================================# DETECCIÓN Y CARGA/SIMULACIÓN DE DATOS# ===============================================# Intentar cargar datos realesdata_path <-paste0(gsub("/docs", "", getwd()), "/data/CONS_C1_2010_22_CLEAN.rds")if(file.exists(data_path)) { data <-readRDS(data_path)} else {# Simular datos con propiedades idénticas a los reales n_sim <-100000# Se puede ajustar este número (original: 223055) data <-simular_datos_chile(n_pacientes = n_sim, seed =2024)}# Función de simplificaciónsimplify_substance_names <-function(x) { x <-as.character(x) x[str_detect(x, "^Sin ")] <-NA x <-str_replace(x, "^Sedantes:.*", "Sedantes") x <-str_replace(x, "^Hipnóticos:.*", "Hipnóticos") x <-str_replace(x, "^Inhalables:.*", "Inhalables") x <-str_replace(x, "^Otros Opioides.*", "Opioides") x <-str_replace(x, "^Otros Estimulantes.*", "Estimulantes") x <-str_replace(x, "^Otros Alucinógenos.*", "Alucinógenos") x <-str_replace(x, "^Éxtasis.*", "Éxtasis/MDMA")return(x)}# Columnas de sustanciascols_sustancias <-c("sustancia_principal", "otras_sustancias_no1", "otras_sustancias_no2", "otras_sustancias_no3")# Aplicar limpiezadata_network <- data %>%select(HASH_KEY, all_of(cols_sustancias)) %>%mutate(across(all_of(cols_sustancias), simplify_substance_names)) %>%filter(!is.na(sustancia_principal))
1.2 INTRODUCCIÓN
Los trastornos por uso de sustancias (TUS) constituyen una de las principales causas de carga de enfermedad y mortalidad evitable a escala global. En 2016 se calculó que más de 100 millones de personas sufrían trastorno por consumo de alcohol y decenas de millones presentaban dependencia de opioides, cannabis o cocaína1. La frecuente comorbilidad psiquiátrica, depresión, trastornos de ansiedad, psicosis o trastornos de personalidad multiplica la severidad clínica y los costes sociosanitarios2. Estudios hospitalarios europeos y norteamericanos muestran que alrededor del 20 % de las admisiones psiquiátricas corresponden a pacientes de sexo femenino con diagnóstico dual, fenómeno que favorece re-ingresos y estancias prolongadas3.
En Chile, las encuestas nacionales sitúan la prevalencia de abuso o dependencia de sustancias entre el 11 % y el 20 %, una de las más elevadas de Latinoamérica. Los registros hospitalarios concuerdan con las cifras internacionales: alrededor de una quinta parte de los internados en psiquiatría presenta un TCS como diagnóstico primario o secundario. Esta convergencia evidencia que la hospitalización psiquiátrica es un desenlace clínico crítico en la trayectoria de las adicciones, razón por la cual identificar sus factores determinantes resulta esencial para planificar intervenciones preventivas, asignar recursos y reducir la carga asistencial2,4–6.
1.3 DESCRIPCIÓN GENERAL DE LOS DATOS
Este es un estudio de cohorte retrospectiva de pacientes adultos en tratamiento por consumo de sustancias, con datos otorgados por el Servicio Nacional para la Prevención y Rehabilitación del Consumo de Drogas y Alcohol de Chile (SENDA) en convenio con el núcleo milenio de ánalisis de políticas públicas de drogas (nDP). La cohorte se construyó vinculando los registros administrativos de los pacientes (n = 223,061 episodios de tratamiento entre 97,698 personas en las 16 regiones del país).
Estos datos incluyen múltiples variables relacionadas al consumo y tratamiento rehabilitador de drogas. Entre estas variables esta la sustancia principal por la cual se trató al paciente y sustancias secundarias (alcohol, pasta base de cocaína, cocaína, marihuana, depresores del SNC u otras sustancias. Tambien está presente el número de reingresos a tratamiento (retratamientos, categorizados en 0, 1, 2, 3 o más reingresos), el tipo de plan de tratamiento (ambulatorio vs. residencial) y el historial clínico de salud mental de los pacientes.
El registro de pacientes en tratamiento se realizó en una plataforma electrónica denominada SISTRAT, que contenía información sociodemográfica, datos sobre el estado de salud y patrones de consumo de sustancias, entre otras variables, además de información sobre el propio tratamiento (p. ej., fecha de ingreso, egreso, tipo de tratamiento). Las base de datos se vincularon de forma determinista mediante un hash de 64 caracteres resultante del cifrado (con un algoritmo SHA-256) del número de identificación único de cada persona.
Código
info_data <-data.frame( Característica =c("Total de registros","Registros con sustancia principal","Período de análisis","Variables de sustancias","Número de sustancias únicas","Promedio de sustancias por paciente","Mediana de sustancias por paciente","Desviación estándar"),Valor =c(format(nrow(data), big.mark =","),format(nrow(data_network), big.mark =","),"2010-2022",as.character(length(cols_sustancias)),as.character(n_distinct(unlist(data_network[cols_sustancias]), na.rm =TRUE)),round(mean(rowSums(!is.na(data_network[cols_sustancias]))), 2),median(rowSums(!is.na(data_network[cols_sustancias]))),round(sd(rowSums(!is.na(data_network[cols_sustancias]))), 2)))info_data %>%kable(format ="html", col.names =c("Característica", "Valor"),align =c("l", "r"),row.names =FALSE) %>%kable_styling(bootstrap_options =c("striped", "hover"),full_width =FALSE)
Características generales de la base de datos
Característica
Valor
Total de registros
223,061
Registros con sustancia principal
223,055
Período de análisis
2010-2022
Variables de sustancias
4
Número de sustancias únicas
22
Promedio de sustancias por paciente
2.2
Mediana de sustancias por paciente
2
Desviación estándar
0.99
2 RED DE CO-OCURRENCIA
2.1 Construcción de la Red
La red de co-ocurrencia es pertinente para identificar patrones de policonsumo y asociaciones entre sustancias que tienden a ser usadas conjuntamente. Esta representación revela combinaciones frecuentes y comunidades de sustancias; en particular, ayuda a localizar “nodos puente” que conectan módulos y que suelen ser dianas clínicas de alto impacto para el cribado y la intervención integrada7–9. Además, varios emparejamientos visibles en la red coinciden con riesgos clínicos conocidos: el uso concurrente de opioides y benzodiacepinas se asocia con un aumento de 2–5 veces del riesgo de sobredosis (mayor en los primeros 90 días de co‑exposición), y más del 90% de las muertes con benzodiacepinas co‑involucran opioides10–12; la co‑ingesta de alcohol y cocaína genera cocaetileno, metabolito más cardiotóxico, por lo que constituye una combinación especialmente peligrosa13; y la co‑involucración de estimulantes con opioides sintéticos (p. ej., fentanilo) ha aumentado de forma marcada en los últimos años, subrayando la necesidad de estrategias específicas para policonsumo14. Metodológicamente, la transformación logarítmica de pesos resulta adecuada en redes con distribuciones de cola pesada, permitiendo visualizar de forma equilibrada tanto las asociaciones muy frecuentes como las moderadas15. Por su parte, filtrar el 60% superior de co‑ocurrencias es coherente con los enfoques de “backbone” en redes ponderadas, que eliminan enlaces débiles para resaltar la estructura multiescala estadísticamente significativa y focalizar los patrones robustos16.
Código
# Crear matriz de co-ocurrenciacreate_cooccurrence_matrix <-function(df) { all_substances <-unique(unlist(df[cols_sustancias])) all_substances <- all_substances[!is.na(all_substances)] all_substances <-sort(all_substances) n_sust <-length(all_substances) co_matrix <-matrix(0, nrow = n_sust, ncol = n_sust,dimnames =list(all_substances, all_substances))for (i in1:nrow(df)) { row_substances <-unlist(df[i, cols_sustancias]) row_substances <- row_substances[!is.na(row_substances)]if (length(row_substances) >1) {for (j in1:(length(row_substances)-1)) {for (k in (j+1):length(row_substances)) { sust1 <- row_substances[j] sust2 <- row_substances[k] co_matrix[sust1, sust2] <- co_matrix[sust1, sust2] +1 co_matrix[sust2, sust1] <- co_matrix[sust2, sust1] +1 } } } }return(co_matrix)}co_matrix <-create_cooccurrence_matrix(data_network)# IMPORTANTE: Solo el 60% de las co-ocurrencias más frecuentes# Obtener todos los valores únicos de co-ocurrencia (excluyendo diagonal)co_values <- co_matrix[upper.tri(co_matrix)]co_values <- co_values[co_values >0]# Calcular el percentil 40 (para mantener el 60% superior)threshold_value <-quantile(co_values, probs =0.40)# Aplicar umbral a la matrizco_matrix_filtered <- co_matrixco_matrix_filtered[co_matrix < threshold_value] <-0# Crear grafo con pesos logarítmicos (ESCALA LOGARÍTMICA)g_full_co <-graph_from_adjacency_matrix( co_matrix_filtered,mode ="undirected",weighted =TRUE,diag =FALSE)# Aplicar transformación logarítmica a los pesosE(g_full_co)$weight_original <-E(g_full_co)$weightE(g_full_co)$weight <-log1p(E(g_full_co)$weight) # ESCALA LOGARÍTMICA: log(1 + weight)# Eliminar vértices aisladosg_co <-delete_vertices(g_full_co, degree(g_full_co) ==0)# Detectar comunidadescommunities_co <-cluster_louvain(g_co)V(g_co)$community <-membership(communities_co)# Información sobre el filtradon_edges_original <-sum(co_matrix >0) /2n_edges_filtered <-sum(co_matrix_filtered >0) /2percent_retained <-round(n_edges_filtered / n_edges_original *100, 1)
Red interactiva de co-ocurrencia (zoom y arrastre habilitados)
3 RED BIPARTITA PACIENTE-SUSTANCIA
3.1 Construcción de la Red
La red bipartita es útil y que representa simultáneamente a los pacientes y a las sustancias, preservando la heterogeneidad individual y, a la vez, permitiendo ver la prevalencia relativa de cada sustancia a través del grado de sus nodos17,18. Al mantener la estructura de dos modos, se evita la pérdida de información y las distorsiones que introducen las proyecciones unipartitas y se pueden aplicar métricas específicas de bipartitos para capturar patrones de alto riesgo19. Además, medidas como la betweenness centrality identifican sustancias que actúan como “puentes” entre grupos de consumidores, clave para detectar rutas de transición en el policonsumo20. Para descubrir perfiles y comunidades de policonsumo se dispone de modularidad bipartita y variantes para redes ponderadas, validadas en la literatura metodológica21,22. La utilidad clínica de los bipartitos está documentada en contextos de salud: redes paciente‑prescriptor y paciente‑hospital han permitido detectar conductas de búsqueda de fármacos y diferenciar patrones según sustancia, evidenciando nodos (sustancias o instituciones) estratégicos para intervención19,Kim2024?. Finalmente, para equilibrar legibilidad y representatividad de la visualización, el muestreo aleatorio simple de 2.500 pacientes facilita la visulización sin introducir sesgos y, usado con fines gráficos, preserva propiedades globales (distribuciones y relaciones clave) de forma aceptable para el análisis exploratorio23,24.
Código
# MUESTREO ALEATORIO SIMPLE DE 2,500 SUJETOS DEL TOTAL# No estratificado por sustancia - cada paciente tiene igual probabilidad de ser seleccionadoset.seed(42)n_sample <-10000sample_patients <-sample(unique(data_network$HASH_KEY), min(n_sample, length(unique(data_network$HASH_KEY))))data_network_sampled <- data_network %>%filter(HASH_KEY %in% sample_patients)# Crear edge list paciente-sustancia con datos muestreadoscreate_bipartite_edgelist <-function(data) { edgelist <-data.frame()for (col in cols_sustancias) { temp <- data %>%select(HASH_KEY, sustancia =!!sym(col)) %>%filter(!is.na(sustancia)) %>%mutate(weight =ifelse(col =="sustancia_principal", 2, 1)) edgelist <-rbind(edgelist, temp) } edgelist <- edgelist %>%group_by(HASH_KEY, sustancia) %>%summarise(weight =sum(weight), .groups ='drop')return(edgelist)}bipartite_edges <-create_bipartite_edgelist(data_network_sampled)# Crear grafo bipartitog_bipartite <-graph_from_data_frame(bipartite_edges, directed =FALSE)V(g_bipartite)$type <-V(g_bipartite)$name %in% bipartite_edges$sustanciaV(g_bipartite)$node_type <-ifelse(V(g_bipartite)$type, "Sustancia", "Paciente")# Detectar comunidadescommunities_bi <-cluster_louvain(g_bipartite)V(g_bipartite)$community <-membership(communities_bi)# Estadísticas básicasn_pacientes <-sum(!V(g_bipartite)$type)n_sustancias <-sum(V(g_bipartite)$type)n_conexiones <-ecount(g_bipartite)
3.2 Métricas de la Red Bipartita
Código
# Calcular grado por tipodegree_all_bi <-degree(g_bipartite)degree_pacientes <- degree_all_bi[!V(g_bipartite)$type]degree_sustancias <- degree_all_bi[V(g_bipartite)$type]# Propiedades estructuralesnetwork_props_bi <-data.frame(Propiedad =c("Número de pacientes (muestra)", "Número de sustancias","Total de conexiones", "Densidad","Conexiones promedio por paciente","Conexiones promedio por sustancia","Componentes conectados", "Modularidad"),Valor =c(format(n_pacientes, big.mark =","), n_sustancias,format(n_conexiones, big.mark =","),round(n_conexiones / (n_pacientes * n_sustancias), 6),round(mean(degree_pacientes), 2),round(mean(degree_sustancias), 0),components(g_bipartite)$no,round(modularity(communities_bi), 3)))network_props_bi %>%kable(format ="html", align =c("l", "r")) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width =FALSE)
Métricas principales de la red bipartita (muestra aleatoria simple de 10000 pacientes)
Propiedad
Valor
Número de pacientes (muestra)
10,000
Número de sustancias
21
Total de conexiones
23,735
Densidad
0.113024
Conexiones promedio por paciente
2.37
Conexiones promedio por sustancia
1130
Componentes conectados
1
Modularidad
0.341
3.3 Visualización Estática Bipartita
Código
# Cargar ggrepelrequire(ggrepel)# Usar el grafo muestreadog_bi_viz <- g_bipartite# Calcular grados para TODOS los nodosdegree_bi_viz <-degree(g_bi_viz)# Calcular tamaños proporcionales al gradonode_sizes_bi <-numeric(vcount(g_bi_viz))# Para sustancias: escala más grandesubstance_nodes <-which(V(g_bi_viz)$type)if(length(substance_nodes) >0) { substance_degrees <- degree_bi_viz[substance_nodes] node_sizes_bi[substance_nodes] <-15+ (substance_degrees /max(substance_degrees)) *30}# Para pacientes: escala más pequeñapatient_nodes <-which(!V(g_bi_viz)$type)if(length(patient_nodes) >0) { patient_degrees <- degree_bi_viz[patient_nodes] node_sizes_bi[patient_nodes] <-0.5+ (patient_degrees /max(patient_degrees)) *5}# ===== OPCIÓN 1: LAYOUT CIRCULAR/RADIAL =====create_radial_layout <-function() { layout_matrix <-matrix(0, vcount(g_bi_viz), 2)# Pacientes en círculo interior (radio menor)if(length(patient_nodes) >0) { angles_patients <-seq(0, 2*pi, length.out =length(patient_nodes) +1)[1:length(patient_nodes)] radius_patients <-3# Radio interior# Variar ligeramente el radio según el grado para crear textura radius_variation <-0.3* (degree_bi_viz[patient_nodes] /max(degree_bi_viz[patient_nodes])) layout_matrix[patient_nodes, 1] <- (radius_patients + radius_variation) *cos(angles_patients) layout_matrix[patient_nodes, 2] <- (radius_patients + radius_variation) *sin(angles_patients) }# Sustancias en círculo exterior (radio mayor)if(length(substance_nodes) >0) {# Ordenar por grado para mejor distribución substance_order <- substance_nodes[order(degree_bi_viz[substance_nodes], decreasing =TRUE)]# Distribuir ángulos con peso según el grado weights <- degree_bi_viz[substance_order] +10 cumulative_weights <-cumsum(weights) angles_substances <- (cumulative_weights /sum(weights)) *2* pi radius_substances <-5# Radio exteriorfor(i inseq_along(substance_order)) { layout_matrix[substance_order[i], 1] <- radius_substances *cos(angles_substances[i]) layout_matrix[substance_order[i], 2] <- radius_substances *sin(angles_substances[i]) } }return(layout_matrix)}# ===== OPCIÓN 2: LAYOUT ARCO =====create_arc_layout <-function() { layout_matrix <-matrix(0, vcount(g_bi_viz), 2)# Pacientes en arco inferiorif(length(patient_nodes) >0) {# Crear un arco de -120 a 120 grados angles <-seq(-2*pi/3, 2*pi/3, length.out =length(patient_nodes)) radius <-4# Agregar variación en radio según grado radius_var <- radius +0.5* (degree_bi_viz[patient_nodes] /max(degree_bi_viz[patient_nodes])) layout_matrix[patient_nodes, 1] <- radius_var *cos(angles) layout_matrix[patient_nodes, 2] <- radius_var *sin(angles) -2# Desplazar hacia abajo }# Sustancias en arco superiorif(length(substance_nodes) >0) { substance_order <- substance_nodes[order(degree_bi_viz[substance_nodes], decreasing =TRUE)]# Distribuir con peso en un arco superior más pequeño weights <-sqrt(degree_bi_viz[substance_order]) +1 positions <-cumsum(c(0, weights[-length(weights)])) positions <- (positions /max(positions)) * pi - pi/2# Arco de -90 a 90 grados radius <-6for(i inseq_along(substance_order)) { layout_matrix[substance_order[i], 1] <- radius *cos(positions[i]) layout_matrix[substance_order[i], 2] <- radius *sin(positions[i]) +3# Desplazar hacia arriba } }return(layout_matrix)}# ===== OPCIÓN 3: LAYOUT GRID EXPANDIDO =====create_expanded_grid_layout <-function() { layout_matrix <-matrix(0, vcount(g_bi_viz), 2)# Pacientes en grid expandido con jitterif(length(patient_nodes) >0) { n_patients <-length(patient_nodes)# Crear una cuadrícula más dispersa n_cols <-ceiling(sqrt(n_patients *2)) # Más columnas para mayor dispersión n_rows <-ceiling(n_patients / n_cols)# Ordenar por grado para agrupar similares patient_order <- patient_nodes[order(degree_bi_viz[patient_nodes], decreasing =TRUE)] idx <-1for(i in1:n_rows) {for(j in1:n_cols) {if(idx <= n_patients) {# Agregar jitter proporcional al grado jitter_strength <-0.2* (degree_bi_viz[patient_order[idx]] /max(degree_bi_viz[patient_nodes])) layout_matrix[patient_order[idx], 1] <- (j - n_cols/2) *0.3+runif(1, -jitter_strength, jitter_strength) layout_matrix[patient_order[idx], 2] <- (i - n_rows/2) *0.3-2# Abajo idx <- idx +1 } } } }# Sustancias en arco superior expandidoif(length(substance_nodes) >0) { substance_order <- substance_nodes[order(degree_bi_viz[substance_nodes], decreasing =TRUE)]# Distribuir en dos filas escalonadas n_sust <-length(substance_order) half <-ceiling(n_sust /2)# Primera fila (sustancias más importantes)for(i in1:min(half, n_sust)) { x_pos <- (i - half/2-0.5) *1.2 layout_matrix[substance_order[i], 1] <- x_pos layout_matrix[substance_order[i], 2] <-3 }# Segunda fila (si hay suficientes sustancias)if(n_sust > half) {for(i in (half+1):n_sust) { x_pos <- ((i-half) - (n_sust-half)/2-0.5) *1.2 layout_matrix[substance_order[i], 1] <- x_pos layout_matrix[substance_order[i], 2] <-2.2 } } }return(layout_matrix)}# SELECCIONAR EL LAYOUT (cambiar según preferencia)# Opción 1: Radial/Circular# layout_matrix <- create_radial_layout()# Opción 2: Arco# layout_matrix <- create_arc_layout()# Opción 3: Grid Expandido (por defecto)set.seed(42) # Para reproducibilidad del jitterlayout_matrix <-create_expanded_grid_layout()# Crear el grafo tidyg_bi_tidy <-as_tbl_graph(g_bi_viz) |>activate(nodes) |>mutate(node_type =ifelse(type, "Sustancia", "Paciente"),node_size_plot = node_sizes_bi,node_degree = degree_bi_viz,label =ifelse(type, name, NA_character_),node_color =ifelse(type, "#e74c3c", "#3498db") )# Crear la visualizaciónp <-ggraph(g_bi_tidy, layout = layout_matrix) +# Edges con curvatura para mejor visualizacióngeom_edge_link(aes(alpha = weight),color ="gray75", width =0.1,show.legend =FALSE) +# Nodosgeom_node_point(aes(size = node_size_plot, color = node_type,alpha = node_type),show.legend =TRUE) +# Etiquetas con ggrepel ggrepel::geom_text_repel(data =function(x) filter(x, !is.na(label)),aes(x = x, y = y, label = label),size =3.5,box.padding =0.5,point.padding =0.3,segment.color ="gray50",segment.size =0.3,segment.alpha =0.5,max.overlaps =50,fontface ="bold",force =2, # Mayor fuerza de repulsiónforce_pull =0.5,max.time =2,max.iter =10000) +scale_size_identity() +scale_color_manual(values =c("Paciente"="#3498db", "Sustancia"="#e74c3c"),name ="Tipo de Nodo") +scale_alpha_manual(values =c("Paciente"=0.5, "Sustancia"=0.9),guide ="none") +scale_edge_alpha_continuous(range =c(0.02, 0.15)) +# Ajustar límites para dar más espaciocoord_cartesian(xlim =c(-8, 8), ylim =c(-5, 6)) +labs(title ="Red Bipartita Paciente-Sustancia",subtitle =paste0("Muestra: ",format(n_pacientes, big.mark =",")," pacientes y ", n_sustancias, " sustancias | ",format(n_conexiones, big.mark =","), " conexiones"),caption ="El tamaño de cada nodo es proporcional a su número de conexiones") +theme_void() +theme(plot.title =element_text(size =18, face ="bold", hjust =0.5),plot.subtitle =element_text(size =13, hjust =0.5, color ="gray40"),plot.caption =element_text(size =10, hjust =0.5, color ="gray50", margin =margin(t =10)),plot.background =element_rect(fill ="white", color =NA),legend.position ="bottom",legend.title =element_text(size =11),legend.text =element_text(size =10))print(p)
Red bipartita paciente-sustancia (muestra aleatoria simple de 2,500 pacientes)
La proyección de la red paciente–sustancia sobre el modo “sustancias” convierte las relaciones indirectas en enlaces directos entre sustancias ponderados por el número de pacientes que comparten, lo que permite analizar la estructura de co‑ocurrencia desde una perspectiva agregada y detectar pares de sustancias que aparecen sistemáticamente en los mismos perfiles de consumo25–27. A diferencia de una red de co‑ocurrencia construida solo con episodios simultáneos, esta proyección integra la experiencia compartida de los pacientes a lo largo del período observado, de modo que las asociaciones reflejan patrones de consumo comunes aunque no coincidan en el tiempo, tal como se justifica en el marco de las redes temporales y en aplicaciones de proyección a partir de datos longitudinales28.
Red interactiva de proyección de sustancias (datos completos)
4.4 Estadísticas Descriptivas de la Red Bipartita Proyectada
Código
# Usar la proyección de la red bipartita como red principalg_analysis <- g_substances# Calcular todas las métricas requeridasmetrics_complete <-data.frame( Métrica =c("Número de nodos", "Número de enlaces","Densidad","Diámetro", "Longitud de camino medio","Clustering global","Clustering medio","Grado medio","Grado máximo","Grado mínimo","Asortatividad","Modularidad"),Valor =c(vcount(g_analysis),ecount(g_analysis),round(edge_density(g_analysis), 4),diameter(g_analysis, weights =NA),round(mean_distance(g_analysis, weights =NA), 3),round(transitivity(g_analysis, type ="global"), 3),round(transitivity(g_analysis, type ="average"), 3),round(mean(degree(g_analysis)), 2),max(degree(g_analysis)),min(degree(g_analysis)),round(assortativity_degree(g_analysis), 3),round(modularity(communities_proj), 3)))metrics_complete %>%kable(format ="html", align =c("l", "r")) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width =FALSE)
Estadísticas descriptivas completas de la red bipartita proyectada
Métrica
Valor
Número de nodos
22.0000
Número de enlaces
188.0000
Densidad
0.8139
Diámetro
2.0000
Longitud de camino medio
1.1860
Clustering global
0.8840
Clustering medio
0.9070
Grado medio
17.0900
Grado máximo
21.0000
Grado mínimo
6.0000
Asortatividad
-0.2600
Modularidad
0.0060
4.5 Medidas de Centralidad Completas (Todas las Sustancias)
Resumen estadístico de las medidas de centralidad para toda la red
Medida
Media
Mediana
Desv. Estándar
Mínimo
Máximo
Coef. Variación
Grado
17.0900
18.0000
4.1500
6.0000
21.0000
0.24
Fuerza
21279.1800
1389.5000
42414.9900
20.0000
130582.0000
1.99
Intermediación
0.0711
0.0271
0.1048
0.0000
0.4198
1.47
Cercanía
0.2566
0.2577
0.0657
0.1243
0.3684
0.26
Eigenvector
0.1735
0.0100
0.3447
0.0001
1.0000
1.99
PageRank
4.5450
1.0320
7.5910
0.6860
24.2730
1.67
4.7 Visualización de la Red según Intermediación
Código
# Visualizar red con tamaño proporcional a intermediacióng_viz <-as_tbl_graph(g_analysis)betweenness_vals <-betweenness(g_analysis, normalized =TRUE)node_sizes_betw <-sqrt(betweenness_vals) *30+5set.seed(42)ggraph(g_viz, layout ='fr') +geom_edge_link(aes(width = weight, alpha = weight),color ="gray50",show.legend =FALSE) +geom_node_point(aes(size = betweenness_vals),color ="#e74c3c",alpha =0.8) +geom_node_text(aes(label = name),size =3,repel =TRUE,force =2) +scale_size_continuous(range =c(2, 15),name ="Intermediación") +scale_edge_width_continuous(range =c(0.2, 2)) +scale_edge_alpha_continuous(range =c(0.2, 0.6)) +labs(title ="Red con Intermediación como Centralidad",subtitle ="Tamaño del nodo proporcional a la intermediación normalizada") +theme_void() +theme(plot.title =element_text(size =14, face ="bold"),plot.subtitle =element_text(size =11))
Visualización de la red según intermediación (betweenness)
5 Pregunta Comparativa
¿La red bipartita proyectada de co-ocurrencia de sustancias presenta patrones de formación aleatorios o sigue principios de organización específicos?
5.1 Hipótesis
H0 (Modelo Aleatorio - Erdős-Rényi): La red se forma por asociaciones aleatorias entre sustancias
H1 (Modelo Scale-Free - Barabási-Albert): La red sigue un proceso de adhesión preferencial con distribución de grado de ley de potencias
5.2 Implementación de Modelos Nulos y Bootstrapping
Código
# Parámetros de la red realn_nodes <-vcount(g_analysis)n_edges <-ecount(g_analysis)avg_degree <-mean(degree(g_analysis))p_connect <-edge_density(g_analysis)# Número de simulaciones para bootstrappingn_sims <-5000# Función para calcular métricascalculate_metrics <-function(g) {c(clustering =transitivity(g, type ="global"),path_length =mean_distance(g, directed =FALSE),diameter =diameter(g, directed =FALSE),assortativity =assortativity_degree(g, directed =FALSE),max_degree =max(degree(g)),degree_variance =var(degree(g)) )}# Métricas de la red realreal_metrics <-calculate_metrics(g_analysis)# Simulaciones Bootstrapset.seed(42)# 1. Modelo Erdős-Rényi (ER)er_metrics <-t(replicate(n_sims, { g_er <-erdos.renyi.game(n_nodes, p_connect, type ="gnp")calculate_metrics(g_er)}))# 2. Modelo Barabási-Albert (BA)ba_metrics <-t(replicate(n_sims, { m_ba <-max(1, round(avg_degree/2)) g_ba <-barabasi.game(n_nodes, m = m_ba, directed =FALSE)calculate_metrics(g_ba)}))# Calcular intervalos de confianza (95%)calculate_ci <-function(metrics_matrix, confidence =0.95) { alpha <-1- confidence lower <-apply(metrics_matrix, 2, quantile, probs = alpha/2, na.rm =TRUE) upper <-apply(metrics_matrix, 2, quantile, probs =1- alpha/2, na.rm =TRUE) mean_val <-colMeans(metrics_matrix, na.rm =TRUE)return(list(mean = mean_val, lower = lower, upper = upper))}er_ci <-calculate_ci(er_metrics)ba_ci <-calculate_ci(ba_metrics)# Crear dataframe con resultados incluyendo intervalos de confianzaresults_df <-data.frame( Métrica =c("Clustering", "Path Length", "Diameter", "Assortativity", "Max Degree", "Degree Variance"),Real = real_metrics,ER_mean = er_ci$mean,ER_CI_lower = er_ci$lower,ER_CI_upper = er_ci$upper,BA_mean = ba_ci$mean,BA_CI_lower = ba_ci$lower,BA_CI_upper = ba_ci$upper)results_df %>%mutate(across(where(is.numeric), round, 3)) %>%mutate(ER =paste0(round(ER_mean, 3), " [", round(ER_CI_lower, 3), ", ", round(ER_CI_upper, 3), "]"),BA =paste0(round(BA_mean, 3), " [", round(BA_CI_lower, 3), ", ", round(BA_CI_upper, 3), "]") ) %>%select(Métrica, Real, ER, BA) %>%kable(format ="html",caption ="Comparación de métricas: Red real vs Modelos nulos (1000 simulaciones, IC 95%)",col.names =c("Métrica", "Red Real", "Erdős-Rényi [IC 95%]", "Barabási-Albert [IC 95%]"),row.names = F) %>%kable_styling(bootstrap_options =c("striped", "hover"), full_width =FALSE)
Comparación de métricas: Red real vs Modelos nulos (1000 simulaciones, IC 95%)
Métrica
Red Real
Erdős-Rényi [IC 95%]
Barabási-Albert [IC 95%]
Clustering
0.884
0.812 [0.757, 0.862]
0.708 [0.694, 0.722]
Path Length
4.208
1.187 [1.139, 1.238]
1.338 [1.338, 1.338]
Diameter
10.000
2 [2, 2]
2 [2, 2]
Assortativity
-0.260
-0.092 [-0.152, -0.029]
-0.13 [-0.224, -0.032]
Max Degree
21.000
20.034 [19, 21]
19.383 [18, 21]
Degree Variance
17.229
3.017 [1.515, 5.16]
11.447 [8.656, 14.563]
Código
# Preparar datos para visualizaciónmetrics_long <-data.frame(Clustering =c(er_metrics[,1], ba_metrics[,1]),PathLength =c(er_metrics[,2], ba_metrics[,2]),Assortativity =c(er_metrics[,4], ba_metrics[,4]),MaxDegree =c(er_metrics[,5], ba_metrics[,5]),Model =rep(c("ER", "BA"), each = n_sims))# Gráficos de distribución con intervalos de confianza como líneas verticalesp1 <-ggplot(metrics_long, aes(x = Clustering, fill = Model)) +geom_density(alpha =0.5) +# Línea vertical para el valor realgeom_vline(xintercept = real_metrics[1], color ="red", linetype ="dashed", size =1.2) +# Intervalos de confianza como líneas verticales para ERgeom_vline(xintercept = er_ci$lower[1], color ="#3498db", linetype ="dotted", size =1) +geom_vline(xintercept = er_ci$upper[1], color ="#3498db", linetype ="dotted", size =1) +# Intervalos de confianza como líneas verticales para BAgeom_vline(xintercept = ba_ci$lower[1], color ="#9b59b6", linetype ="dotted", size =1) +geom_vline(xintercept = ba_ci$upper[1], color ="#9b59b6", linetype ="dotted", size =1) +scale_fill_manual(values =c("#3498db", "#9b59b6")) +labs(title ="Clustering", x ="Valor", y ="Densidad") +theme_minimal() +annotate("text", x = real_metrics[1], y =0, label ="Real", color ="red", vjust =-1)p2 <-ggplot(metrics_long, aes(x = PathLength, fill = Model)) +geom_density(alpha =0.5) +geom_vline(xintercept = real_metrics[2], color ="red", linetype ="dashed", size =1.2) +# Intervalos de confianza como líneas verticales para ERgeom_vline(xintercept = er_ci$lower[2], color ="#3498db", linetype ="dotted", size =1) +geom_vline(xintercept = er_ci$upper[2], color ="#3498db", linetype ="dotted", size =1) +# Intervalos de confianza como líneas verticales para BAgeom_vline(xintercept = ba_ci$lower[2], color ="#9b59b6", linetype ="dotted", size =1) +geom_vline(xintercept = ba_ci$upper[2], color ="#9b59b6", linetype ="dotted", size =1) +scale_fill_manual(values =c("#3498db", "#9b59b6")) +labs(title ="Path Length", x ="Valor", y ="Densidad") +theme_minimal() +annotate("text", x = real_metrics[2], y =0, label ="Real", color ="red", vjust =-1)p3 <-ggplot(metrics_long, aes(x = Assortativity, fill = Model)) +geom_density(alpha =0.5) +geom_vline(xintercept = real_metrics[4], color ="red", linetype ="dashed", size =1.2) +# Intervalos de confianza como líneas verticales para ERgeom_vline(xintercept = er_ci$lower[4], color ="#3498db", linetype ="dotted", size =1) +geom_vline(xintercept = er_ci$upper[4], color ="#3498db", linetype ="dotted", size =1) +# Intervalos de confianza como líneas verticales para BAgeom_vline(xintercept = ba_ci$lower[4], color ="#9b59b6", linetype ="dotted", size =1) +geom_vline(xintercept = ba_ci$upper[4], color ="#9b59b6", linetype ="dotted", size =1) +scale_fill_manual(values =c("#3498db", "#9b59b6")) +labs(title ="Assortativity", x ="Valor", y ="Densidad") +theme_minimal() +annotate("text", x = real_metrics[4], y =0, label ="Real", color ="red", vjust =-1)p4 <-ggplot(metrics_long, aes(x = MaxDegree, fill = Model)) +geom_density(alpha =0.5) +geom_vline(xintercept = real_metrics[5], color ="red", linetype ="dashed", size =1.2) +# Intervalos de confianza como líneas verticales para ERgeom_vline(xintercept = er_ci$lower[5], color ="#3498db", linetype ="dotted", size =1) +geom_vline(xintercept = er_ci$upper[5], color ="#3498db", linetype ="dotted", size =1) +# Intervalos de confianza como líneas verticales para BAgeom_vline(xintercept = ba_ci$lower[5], color ="#9b59b6", linetype ="dotted", size =1) +geom_vline(xintercept = ba_ci$upper[5], color ="#9b59b6", linetype ="dotted", size =1) +scale_fill_manual(values =c("#3498db", "#9b59b6")) +labs(title ="Max Degree", x ="Valor", y ="Densidad") +theme_minimal() +annotate("text", x = real_metrics[5], y =0, label ="Real", color ="red", vjust =-1)# Combinar los gráficos(p1 + p2) / (p3 + p4) +plot_layout(guides ="collect") &theme(legend.position ="bottom")
Distribuciones bootstrap de métricas clave con intervalos de confianza
Código
# Generar distribuciones de grado con IC para datos realesset.seed(42)n_bootstrap <-5000# Función para calcular distribución de gradocalc_degree_dist <-function(g) { deg <-degree(g) deg_table <-table(deg)data.frame(k =as.numeric(names(deg_table)),p_k =as.numeric(deg_table) /sum(deg_table) )}# Bootstrap para la red real - resampleando nodosreal_dist_bootstrap <-list()nodes_total <-vcount(g_analysis)for(i in1:n_bootstrap) {# Resamplear nodos con reemplazo sampled_nodes <-sample(1:nodes_total, nodes_total, replace =TRUE) g_sampled <-induced_subgraph(g_analysis, sampled_nodes)# Calcular distribución del grafo resampleadoif(vcount(g_sampled) >0&&ecount(g_sampled) >0) { real_dist_bootstrap[[i]] <-calc_degree_dist(g_sampled) }}# Calcular distribución promedio e IC para la red realaggregate_real_distribution <-function(dist_list) { all_k <-sort(unique(unlist(lapply(dist_list, function(x) x$k))))# Crear matriz para almacenar todas las p_k n_valid <-length(dist_list) p_k_matrix <-matrix(0, nrow = n_valid, ncol =length(all_k))colnames(p_k_matrix) <-as.character(all_k)for(i in1:n_valid) { dist <- dist_list[[i]]for(j in1:nrow(dist)) { col_idx <-which(all_k == dist$k[j]) p_k_matrix[i, col_idx] <- dist$p_k[j] } }# Calcular media y desviación estándar mean_p_k <-colMeans(p_k_matrix, na.rm =TRUE) sd_p_k <-apply(p_k_matrix, 2, sd, na.rm =TRUE)# IC usando ±1 desviación estándar (como en el ejemplo del paper)data.frame(k = all_k,p_k = mean_p_k,p_k_lower =pmax(mean_p_k - sd_p_k, 0), # No puede ser negativop_k_upper = mean_p_k + sd_p_k,Model ="Real" )}# Calcular IC para datos realesdd_real_ci <-aggregate_real_distribution(real_dist_bootstrap)# Generar distribuciones teóricas para ER y BA (sin IC, solo línea teórica)# Para ER: usar la fórmula teórica de distribución binomialn_nodes <-vcount(g_analysis)p_connect <-edge_density(g_analysis)k_range <-0:50# Distribución teórica ER (binomial)dd_er_theory <-data.frame(k = k_range,p_k =dbinom(k_range, n_nodes -1, p_connect),p_k_lower =NA,p_k_upper =NA,Model ="ER")# Para BA: generar una muestra grande para obtener la distribución esperadam_ba <-max(1, round(mean(degree(g_analysis))/2))g_ba_large <-barabasi.game(n_nodes *10, m = m_ba, directed =FALSE)dd_ba_temp <-calc_degree_dist(g_ba_large)# Interpolar para tener valores en el mismo rango de kdd_ba_theory <-data.frame(k = k_range,p_k =approx(dd_ba_temp$k, dd_ba_temp$p_k, xout = k_range, rule =2)$y,p_k_lower =NA,p_k_upper =NA,Model ="BA")dd_ba_theory$p_k[is.na(dd_ba_theory$p_k)] <-0# Combinar todos los datosdd_all <-rbind(dd_real_ci, dd_er_theory, dd_ba_theory)# Escala lineal con IC solo para datos realesp1 <-ggplot(dd_all, aes(x = k, y = p_k, color = Model)) +# Banda de confianza solo para datos realesgeom_ribbon(data = dd_real_ci,aes(ymin = p_k_lower, ymax = p_k_upper, fill = Model), alpha =0.3, color =NA) +# Líneas para todos los modelosgeom_line(size =1.2, alpha =0.9) +scale_color_manual(values =c("Real"="#e74c3c", "ER"="#2ecc71", "BA"="#9b59b6"),labels =c("Real"="Datos reales (media)", "ER"="ER (teórico)", "BA"="BA (teórico)")) +scale_fill_manual(values =c("Real"="#e74c3c"),labels =c("Real"="Datos reales (±1 sd)"),guide =guide_legend(override.aes =list(alpha =0.3))) +labs(title ="(a) Escala lineal",x ="k", y =expression(p[k])) +theme_minimal() +theme(plot.title =element_text(size =14, face ="bold"),legend.position ="bottom",legend.title =element_blank()) +coord_cartesian(xlim =c(0, 100))# Escala log-log con IC solo para datos reales# Filtrar valores positivosdd_log <- dd_all %>%filter(k >0, p_k >0)dd_real_log <- dd_real_ci %>%filter(k >0, p_k >0) %>%mutate(p_k_lower_safe =pmax(p_k_lower, 1e-6),p_k_upper_safe = p_k_upper )p2 <-ggplot(dd_log, aes(x = k, y = p_k, color = Model)) +# Banda de confianza solo para datos realesgeom_ribbon(data = dd_real_log,aes(ymin = p_k_lower_safe, ymax = p_k_upper_safe, fill = Model), alpha =0.3, color =NA) +# Líneas para todos los modelosgeom_line(size =1.2, alpha =0.9) +# Agregar línea de tendencia para BA (ley de potencia)geom_smooth(data = dd_log[dd_log$Model =="BA"& dd_log$k >10,],method ="lm", formula = y ~ x, se =FALSE,linetype ="dashed", size =0.8, alpha =0.5) +scale_x_log10(breaks =c(1, 10, 30, 50),labels =c("1", "10", "30", "50")) +scale_y_log10(breaks =c(1e-4, 1e-3, 1e-2, 1e-1, 1),labels = scales::trans_format("log10", scales::math_format(10^.x))) +scale_color_manual(values =c("Real"="#e74c3c", "ER"="#2ecc71", "BA"="#9b59b6"),labels =c("Real"="Datos reales (media)", "ER"="ER (teórico)", "BA"="BA (teórico)")) +scale_fill_manual(values =c("Real"="#e74c3c"),labels =c("Real"="Datos reales (±1 sd)"),guide =guide_legend(override.aes =list(alpha =0.3))) +labs(title ="(b) Escala log-log",x ="k (log)", y =expression(p[k]~"(log)")) +theme_minimal() +theme(plot.title =element_text(size =14, face ="bold"),legend.position ="bottom",legend.title =element_blank(),panel.grid.minor =element_line(size =0.1))# Combinar ambos gráficosp1 + p2 +plot_layout(guides ="collect") &theme(legend.position ="bottom")
Comparación de distribuciones de grado con modelos teóricos e intervalos de confianza
1. Volkow, N. D., & Blanco, C. (2023). Substance use disorders: a comprehensive update of classification, epidemiology, neurobiology, clinical aspects, treatment and prevention. World Psychiatry, 22(2), 203-229. https://doi.org/10.1002/wps.21073
2. Connery, H. S., McHugh, R. K., Reilly, M., Shin, S., & Greenfield, S. F. (2020). Substance Use Disorders in Global Mental Health Delivery: Epidemiology, Treatment Gap, and Implementation of Evidence-Based Treatments. Harvard Review of Psychiatry, 28(5), 316-327. https://doi.org/10.1097/HRP.0000000000000271
3. Gómez-Sánchez-Lafuente, C., Guzman-Parra, J., Suarez-Perez, J., Bordallo-Aragon, A., Rodriguez-de-Fonseca, F., & Mayoral-Cleries, F. (2022). Trends in Psychiatric Hospitalizations of Patients With Dual Diagnosis in Spain. Journal of Dual Diagnosis, 18(2), 92-100. https://doi.org/10.1080/15504263.2022.2053770
4. Saxena, S., Thornicroft, G., Knapp, J., & Whiteford, M. (2007). Resources for mental health: scarcity, inequity, and inefficiency. World Psychiatry, 6(1), 1-10. https://doi.org/10.1002/wps.21073
5. Gómez-Sánchez-Lafuente, C., Guzman-Parra, J., Suarez-Perez, J., Mayoral-Cleries, F., & Rodriguez-de-Fonseca, F. (2016). Dual Diagnosis in Spain: Prevalence, Sociodemographic Profile, and Psychiatric Comorbidity in a Sample of Patients Admitted to Psychiatric Inpatient Units. Journal of Dual Diagnosis, 12(3-4), 249-258. https://doi.org/10.1080/15504263.2016.1220207
6. Rojas, G., Gaete, M., Guajardo, M., Martínez, M., Martínez, M., Fritsch, R., & Araya, R. (2002). Prevalencia de trastornos psiquiátricos en pacientes hospitalizados en un hospital general. Revista Médica de Chile, 130(6), 689-696. https://doi.org/10.4067/S0034-98872002000600008
7. Jones, P. J., Ma, R., & McNally, R. J. (2021). Bridge Centrality: A Network Approach to Understanding Comorbidity. Multivariate Behavioral Research, 56(2), 353-367. https://doi.org/10.1080/00273171.2019.1614898
9. López-Toro, E., Wolf, C. J. H., González, R. A., Brink, W. van den, Schellekens, A., & Vélez-Pastrana, M. C. (2022). Network Analysis of DSM Symptoms of Substance Use Disorders and Frequently Co-Occurring Mental Disorders in Patients with Substance Use Disorder Who Seek Treatment. Journal of Clinical Medicine, 11(10), 2883. https://doi.org/10.3390/jcm11102883
10. Sun, E. C., Dixit, A., Humphreys, K., Darnall, B. D., Baker, L. C., & Mackey, S. (2017). Association between concurrent use of prescription opioids and benzodiazepines and overdose: retrospective analysis. BMJ, 356, j760. https://doi.org/10.1136/bmj.j760
11. Hernandez, I., He, M., Brooks, M. M., & Zhang, Y. (2018). Exposure-Response Association Between Concurrent Opioid and Benzodiazepine Use and Risk of Opioid-Related Overdose in Medicare Part D Beneficiaries. JAMA Network Open, 1(2), e180919. https://doi.org/10.1001/jamanetworkopen.2018.0919
12. Liu, S., O’Donnell, J., Gladden, R. M., McGlone, L., & Chowdhury, F. (2021). Trends in Nonfatal and Fatal Overdoses Involving Benzodiazepines — 38 States and the District of Columbia, 2019–2020. MMWR. Morbidity and Mortality Weekly Report, 70(34), 1136-1141. https://doi.org/10.15585/mmwr.mm7034a2
13. Pennings, E. J. M., Leccese, A. P., & Wolff, F. A. de. (2002). Effects of Concurrent Use of Alcohol and Cocaine. Addiction, 97(7), 773-783. https://doi.org/10.1046/j.1360-0443.2002.00158.x
14. Mattson, C. L., Tanz, L. J., Quinn, K., Kariisa, M., Patel, R., & Davis, N. L. (2021). Trends and Geographic Patterns in Drug and Synthetic Opioid Involvement in Overdose Deaths — United States, 2013–2019. MMWR. Morbidity and Mortality Weekly Report, 70(6), 202-207. https://doi.org/10.15585/mmwr.mm7006a4
15. Clauset, A., Shalizi, C. R., & Newman, M. E. J. (2009). Power-Law Distributions in Empirical Data. SIAM Review, 51(4), 661-703. https://doi.org/10.1137/070710111
16. Serrano, M. Á., Boguña, M., & Vespignani, A. (2009). Extracting the Multiscale Backbone of Complex Weighted Networks. Proceedings of the National Academy of Sciences of the United States of America, 106(16), 6483-6488. https://doi.org/10.1073/pnas.0808904106
17. Wang, T., Bendayan, R., Msosa, Y., Pritchard, M., Roberts, A., Stewart, R., & Dobson, R. (2022). Patient-centric characterization of multimorbidity trajectories in patients with severe mental illnesses: A temporal bipartite network modeling approach. Journal of Biomedical Informatics, 127, 104010. https://doi.org/10.1016/j.jbi.2022.104010
19. Yang, K.-C., Aronson, B., Odabas, M., Ahn, Y.-Y., & Perry, B. L. (2022). Comparing measures of centrality in bipartite patient–prescriber networks: A study of drug seeking for opioid analgesics. PLOS ONE, 17(8), e0273569. https://doi.org/10.1371/journal.pone.0273569
20. Pavlopoulos, G. A., Kontou, P. I., Pavlopoulou, A., Bouyioukos, C., Markou, E., & Bagos, P. G. (2018). Bipartite graphs in systems biology and medicine: a survey of methods and applications. GigaScience, 7(4), giy014. https://doi.org/10.1093/gigascience/giy014
21. Dormann, C. F., & Strauss, R. (2014). A method for detecting modules in quantitative bipartite networks. Methods in Ecology and Evolution, 5(1), 90-98. https://doi.org/10.1111/2041-210X.12139
22. Beckett, S. J. (2016). Improved community detection in weighted bipartite networks. Royal Society Open Science, 3(1), 140536. https://doi.org/10.1098/rsos.140536
23. Wu, Y., Cao, N., Archambault, D., Shen, Q., Qu, H., & Cui, W. (2017). Evaluation of Graph Sampling: A Visualization Perspective. IEEE Transactions on Visualization and Computer Graphics, 23(1), 401-410. https://doi.org/10.1109/TVCG.2016.2598867
25. Latapy, M., Magnien, C., & Del Vecchio, N. (2008). Basic notions for the analysis of large two-mode networks. Social Networks, 30(1), 31-48. https://doi.org/10.1016/j.socnet.2007.04.006
26. Newman, M. E. J. (2001). Scientific collaboration networks. I. Network construction and fundamental results. Physical Review E, 64(1), 016131. https://doi.org/10.1103/PhysRevE.64.016131
27. Guillaume, J.-L., & Latapy, M. (2006). Bipartite graphs as models of complex networks. Physica A: Statistical Mechanics and its Applications, 371(2), 795-813. https://doi.org/10.1016/j.physa.2006.04.047
